unused functions/simulation_sum.R

In BriceonWiley/IntegratedLikelihood.R:

simulation_sum <- function(N, p = 0.1, phi = 0.1, theta = 0.1,
                           k = 0.5, alpha = 0.05, len = 1e4) {
  # browser()
  pip <- p * (1 - theta) + (1 - p) * phi
  n <- 0.1 * N
  M <- 0.9 * N

  if(length(k) > 1){
    results <-
      rerun(len, {
        nprobs <- c((1 - p) * (1 - phi), (1 - p) * phi, p * theta, p * (1 - theta))
        ns <- rmultinom(1, n, nprobs)
        XY <- rmultinom(1, M, c(pip, 1 - pip))
        # mW <- mWald(ns[1], ns[2], ns[3], ns[4], XY[1], XY[2], alpha)
        # mWwidth <- mW[2] - mW[1]
        # mWcov <- p > mW[1] && p < mW[2]
        IL <- intLik(ns[1], ns[2], ns[3], ns[4], XY[1], XY[2], alpha)
        ILwidth <- IL[2] - IL[1]
        ILcov <- p > IL[1] && p < IL[2]
        aIL <- k %>%
          map(~adj_intLik(ns[1], ns[2], ns[3], ns[4], XY[1], XY[2], k = .x, alpha))
        aILwidth <- map(aIL, ~unname(diff(.x)))
        names(aILwidth) <- glue("width_k_{k}")
        aILwidth %<>% unlist() %>% t() %>% as_tibble()
        aILcov <- map(aIL, ~between(p, min(.x), max(.x)))
        names(aILcov) <- glue("cov_k_{k}")
        aILcov %<>%  unlist() %>% t() %>% as_tibble()
        bind_cols(
          # "mWwidth" = mWwidth, "mWcov" = mWcov,
          "ILwidth" = ILwidth, "ILcov" = ILcov,
          aILwidth, aILcov
        )
      }) %>%
      bind_rows() %>%
      summarize_all(mean, na.rm = T)

    # mWres <- c(
    #   "Average Width" = results$mWwidth,
    #   "Coverage Probability" = results$mWcov
    # )
    ILres <- c(
      "Average Width" = results$ILwidth,
      "Coverage Probability" = results$ILcov
    )
    aILres <- matrix(results[, -(1:2)], length(k))
    row.names(aILres) <- glue("adjusted Integrated Likelihood k = {k}")
    colnames(aILres) <- c("Average Width", "Coverage Probability")

    list(
      "counts" = c(
        "N" = N,
        "M" = M,
        "n" = n,
        "B" = len
      ),
      "proportions" = c(
        "p" = p,
        "phi" = phi,
        "theta" = theta,
        "k" = k
      ),
      "results" = rbind(
        # "mWald" = mWres,
        "Integrated Likelihood" = ILres,
        aILres
      )
    )
  } else{
    results <-
      rerun(len, {
        nprobs <- c((1 - p) * (1 - phi), (1 - p) * phi, p * theta, p * (1 - theta))
        ns <- rmultinom(1, n, nprobs)
        XY <- rmultinom(1, M, c(pip, 1 - pip))
        # mW <- mWald(ns[1], ns[2], ns[3], ns[4], XY[1], XY[2], alpha)
        # mWwidth <- mW[2] - mW[1]
        # mWcov <- p > mW[1] && p < mW[2]
        IL <- intLik(ns[1], ns[2], ns[3], ns[4], XY[1], XY[2], alpha)
        ILwidth <- IL[2] - IL[1]
        ILcov <- p > IL[1] && p < IL[2]
        aIL <- adj_intLik(ns[1], ns[2], ns[3], ns[4], XY[1], XY[2], k = k, alpha)
        aILwidth <- aIL[2] - aIL[1]
        aILcov <- p > aIL[1] && p < aIL[2]
        data.frame(
          # "mWwidth" = mWwidth, "mWcov" = mWcov,
          "ILwidth" = ILwidth, "ILcov" = ILcov,
          "aILwidth" = aILwidth, "aILcov" = aILcov
        )
      }) %>%
      bind_rows() %>%
      summarize_all(mean, na.rm = T)

    # mWres <- c(
    #   "Average Width" = results$mWwidth,
    #   "Coverage Probability" = results$mWcov
    # )
    ILres <- c(
      "Average Width" = results$ILwidth,
      "Coverage Probability" = results$ILcov
    )
    aILres <- c(
      "Average Width" = results$aILwidth,
      "Coverage Probability" = results$aILcov
    )

    list(
      "counts" = c(
        "N" = N,
        "M" = M,
        "n" = n,
        "B" = len
      ),
      "proportions" = c(
        "p" = p,
        "phi" = phi,
        "theta" = theta,
        "k" = k
      ),
      "results" = rbind(
        # "mWald" = mWres,
        "Integrated Likelihood" = ILres,
        "adjusted Integrated Likelihood" = aILres
      )
    )
  }
}